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Abstract. 

Using the Ricci and scalar curvatures of the configuration 
manifold of gravitational TV-body systems, we study the ex- 
ponential instability in their trajectories. It is found that the 
exponentiation time-scale for isotropic Plummer spheres varies 
very little with particle number if the softening is small. Large 
softening on the other hand has a marked effect and, if large 
enough, can cause the curvatures to become positive. This last 
result confirms the previous observations for self gravitating 
sheets and suggests that the qualitative behaviour of large- A'^ 
and continuum systems may be different, and that their equiv- 
alence is only obtained in the limit of infinite A*' and finite 
softening. It is also found that the presence of a large fraction 
of the kinetic energy in rotational motion increases the expo- 
nentiation time-scales significantly — an effect that should be 
expected given the regular nature of nearly circular motion. In 
the light of the results of this and of previous studies, it is sug- 
gested that the exponential instability may arise from low order 
resonances between the period of the variation of the gravita- 
tional field due to distant encounters and the orbital period of 
a test particle. For periods long compared to the exponenti- 
ation time but short compared to the diffusion time-scales of 
the action variables, the standard picture of coUisionless dy- 
namics may be valid in an averaged sense — nevertheless this 
time interval need not coincide with that predicted by stan- 
dard relaxation theory. Instead it is suggested that, at least 
for systems with well defined final states, the relaxation time 
should scale as ~ N^^^. 
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1. Introduction 

It is well documented in numerical simulations that A^-body 
gravitational systems display an exponential instability with 
respect to small changes in the initial conditions (e.g.. Miller 
1964; Goodman et al. 1993; Kandrup et al. 1994). This insta- 
bility not only appears when the linearised dynamics is studied 
but also when the full nonlinear evolution of two originally sim- 
ilar systems is considered. It does not appear therefore to be 
simply a product of linearisation of the equations of motion 
— a linearisation instability. It is possible however (El-Zant 



1996a) that the exponential instability observed in short time 
calculations of the divergence of temporal states may result 
from phase mixing. However there is another way in which the 
instability manifests itself, namely through the predominantly 
negative two dimensional curvature of the configuration space 
of A^-body systems (Gurzadyan & Savvidy 1984,1986; Kan- 
drup 1990a, 1990b). This property arises from the qualitative 
phase space structure of a system and cannot be explained 
away so easily. It necessarily implies instability normal and 
not just along the phase space^ trajectory of a system (a prop- 
erty known in dynamical literature as transversality: e.g., Ru- 
elle 1989), which implies qualitatively different behaviour from 
that of regular systems. Since, if real, such an effect could have 
far reaching implications on the evolution of gravitational sys- 
tems, it is natural to enquire as to how the predicted instabil- 
ity correlates with various parameters of A'^-body systems in 
an attempt to uncover its origin and implications. Since not 
all two dimensional curvatures are always negative during the 
evolution of these systems, and in any case, these depend on 
the full Riemann tensor (difficult to calculate for all but the 
smallest systems) an averaged chaos indicator still relying on 
the geometric approach has to be used. 

Such a method, involving the Ricci (or mean) curvature, 
was devised by Gurzadyan & Kocharyan (1987) and was ap- 
plied to numerically integrated A^-body systems by El-Zant 
(1996a) and El-Zant & Gurzadyan (1997). In these papers it 
was shown that the stability of motion of A'^-body gravitational 
systems as described by the aforementioned method correlated 
strongly with various parameters of a gravitational system. 
For although the exponential instability (as evidenced by the 
negative Ricci curvature) exists for most initial conditions, it 
was observed for example that this instability was more pro- 
nounced when clear macroscopic (e.g., collective plasma type) 
instabilities were present and in situations where one expects 
a faster evolution rate (e.g., in the presence of central concen- 
tration which acceUerates the gravothermal evolution of grav- 
itational systems towards still more concentrated states with 
higher thermodynamic entropy: Saslaw 1985). On the other 
hand, the predicted instability was much weaker for regular 
systems (e.g. for ones in uniform rotation before macroscopic 
instability starts to have a significant effect). Moreover, it was 



^The term phase space in this paper will refer to the 6 A'^ 
dimensional phase space unless otherwise indicated. 
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found (El-Zant 1996a) that the negativity of the Ricci curva- 
ture is not a result of contributions duo to close encounters but 
in fact in spite of them — the curvature became always neg- 
ative only when the fluctuations due to the closest encounters 
were removed. In fact, for a system in statistical equilibrium, 
the Ricci curvature is constant (up to random fluctuations and 
neglecting close encounters) and negative and mainly deter- 
mined by the first derivatives of the potential. This suggests 
that the instability arising from this property is unlikely to be 
the result of the linearisation of a system containing fluctuat- 
ing forces. In other words, it is likely to be the manifestation 
of chaotic behaviour characteristic of generic A'^-body gravi- 
tational systems. It is often assumed that if this is the case 
then the cS'ect of the instability mainly concerns the accu- 
racy of numerical simulations of A''-body systems (e.g., Good- 
man et al. 1993; Miller 1994). This problem is also mentioned 
in many papers dealing with numerical simulation techniques 
(e.g., Aarseth 1996; Barnes & Hut 1989). Chaos however does 
not only afltect rmmcrical simulations, it usually comes with 
other deep implications. The most important is that it leads to 
diffusion in the action variables which determine the physical 
state of a system. This diffusion time-scale need not be deter- 
mined by the classical two body theory — even if it is long com- 
pared to the exponentiation time since this theory assumes 
an integrable system which remains so under perturbations 
due to discreteness. On the other hand, systems with negative 
curvature (whether constant or not) usually have strong sta- 
tistical properties, very different from near integrable systems 
(e.g., Pesin 1989). It is therefore very plausible that there may 
be, contrary to what is often assumed, a physical meaning to 
the exponential instability. This meaning is however not yet 
completely clear. 

In this paper we continue our investigation of the stability 
of motion of iV-body systems using the Ricci and scalar curva- 
tures of the configuration manifold of generalised coordinates 
(the method is briefly described in Section 2). We will be inter- 
ested here in the case of spherical systems. The main questions 
we would like to ask are as follows. First, it has been observed 
(both in direct simulations and by using the Ricci curvature 
method) that the exponentiation time-scales are rather small 
(usually a fraction of a dynamical time). If this time-scale does 
not increase with A'^ then one must conclude that it cannot 
be directly related to the diffusion of the action variables in 
a system in virial equilibrium — and therefore any physical 
meaning will have to be more subtle. In Section 3 we exam- 
ine the variation of that time-scale for N up to 15 000. For flat 
sheet like systems it was found (El-Zant 1996a) that significant 
softening destroyed the negative Ricci curvature of the config- 
uration space, meaning that the structure of the phase-space 
should be different and that the motion should be more reg- 
ular. If this is the case, then the exponential instability must 
arise from discreteness effects and there has to be a discontirm- 
ity in the transition between the phase space structure of large 
(but finite) A^-body gravitational systems and their continuum 
counterparts. In Section 4 we study the effect of softening in 
spherical Plummer models and its possible implications in some 
detail. Another problem is concerned with the fact that many 
galaxies are observed to be in fairly regular states where most 
of the kinetic energy is in the form of ordered rotational mo- 
tion. Again, a small exponentiation time-scale (if directly in- 
terpreted) is incompatible with such motion — especially if the 
divergence takes place from most initial conditions, in which 



case there is no room for "stable chaos" or trapping between 
KAM tori (as for example suggested by Gouda et al. 1994). It 
is therefore important to see if the presence of large amount 
of rotational motion does influence the relaxation time-scales 
significantly. This effect is studied in Section 5. In the final 
section we summarise the results and describe a possible inter- 
pretation of the origin and effect of the exponential instability 
of Af-body systems and how this can be tested. 

2. Method 

The study of the stability properties of gravitational A'^-body 
systems using the Ricci curvature was initiated by Gurzadyan 

&: Kocharyan (1987). Details of how the Ricci curvature can 
be used to study the stability of A^-body trajectories are given 
in El-Zant (1996a). We just mention here a few points that are 
essential to the interpretation of the results of the following 
sections. 

Through the Lagrangian formulation of dynamics, it is pos- 
sible to reduce the study of the stability of motion of A^- 
body systems to that of the geometry of the corresponding 
Lagrangian manifold. When this is done, the Gaussian cur- 
vatures on any two dimensional directions normal to the mo- 
tion are found to be mostly (but not always) negative (e.g., 
Gurzadyan & Savvidy 1986; Kandrup 1990a). Since for large 
A' the probability any of the two dimensional curvatures being 
positive becomes exceedingly small, it is possible to replace the 
full set of two dimensional curvatures by their average over the 
3A'^ — 1 directions normal to the geodesic velocity vector u . 
This quantity also happens to coincide with the Ricci curvature 
of the manifold. In Cartesian coordinates in the enveloping 3A'^ 
configuration space this is given by 



ru = 3A 
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with A = ^^^^ . Here W is the kinetic energy calculated as 
a function of the potential energy while VW = ^ Wi and 
V^W — ^ Wii. For a TV-body system, the implied summation 
would be over i,j = 1, SA'^. If we label by a, b and c the particle 
numbers (which run from 1 to A'^) and by k and / the three 
Cartesian coordinates of a particle, then for particles with same 
mass m and with G = 1 we have 
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if a = 6. In these equations 



rib = {ribf + {rlbf + {rtbf + b^, 



(2) 



(3) 



(4) 



(5) 



where 6 is a softening parameter and r*,, — 

A negative Ricci curvature can be interpreted to imply that 
a system will, in general, display exponential instability to ran- 
dom perturbations (El-Zant 1996b). The associated exponen- 
tiation time-scale is then given by 
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where a bar denotes a phase space average (can be either static 
or dynamic). 

When the two dimensional curvatures are averaged over all 
geodesies that pass by a given point as well as over the direc- 
tions normal to the velocity vector at that point, one obtains 
the curvature scalar 

Ji=^A;u,n = ^ru. (7) 

u,n u 

This quantity turns out to be always negative for N > 2. This 
means that a iV-body system will, in general, display expo- 
nential instability when randomly perturbed and for random 
distribution in velocity space. An cissociated instability time- 
scale can also be inferred: One replaces Vu in Eq. (6) by R and 
3N by {3N)^. Clearly the time-scale derived from the scalar 
curvature can be equal to that inferred from the Ricci curvar 
ture only if no ordered motion is present. We will see that this 
is indeed the case. 

It is the exponentiation time-scales that will be of interest 
in the coming sections. We will be averaging them either over 
(pseudo)random realizations of iV-body sytems or over time- 
averages for numerically integrated systems. In the latter case, 
we use the Aarseth (1996) NB0DY2 code which is a direct 
sununation code applying the Ahmad & Cohen (1973) neigh- 
bour scheme and individual time-steps for each of the particles 
in the simulation. These refinements speed up the integration 
considerably, while essentially maintaining the accuracy and 
simplicity of direct summation codes. Since we are still in the 
exploratory stages of applying and testing geometric methods 
to gravitational systems, it is prudent to study the behaviour 
of the curvatures when calculated in an accurate and straight- 
forward manner before integrating it into large-A'^ codes. The 
parameters of the code are fixed at the same values used in 
El-Zant & Gurzadyan (1997). 

Finally, it is important to note here that the negativity of 
the curvatures is only one mechanism by which chaos can occur 
according to the geometric formulation of dynamics. Another 
mechanism is provided by parametric instability as pointed out 
by Pettini (1993) and Cerruti-Sola & Pettini (1995). However, 
when A'' is large and the system considered is near virial equilib- 
rium, the mean potential energy does not undergo significant 
fiuctuations. In that case the second and third terms of Eq. 
(26) of the latter paper, which depend on the time derivatives 
of the potential energy, are then very small. This leaves nega- 
tive curvature as the only effective mechanism for instability. 

3. Behaviour of the curvatures as particle numbers 

change 

In this section we look at the effect of varying the particle imm- 
ber on the exponentiation time-scale of isotropic equilibrium 
Plummer models prepared using the method of Aarseth et al. 
(1974) and scaled to the units of Heggie & Mathieu (1986) by 
keeping the total mass and the gravitational constant equal to 
unity and the total energy fixed at E = —0.25. In this case the 
mean crossing time Tc is equal to 2\/2 time units. 

Direct summation routines are not well suited for integra- 
tion of large AT-body systems (except if used in conjunction 
with special hardware like GRAPE architecture which was not 



available to the author) and NBODY2 becomes very slow for 
particle numbers exceeding a few thousand as to prevent sys- 
tematic examination of the dynamics for greater particle num- 
bers. Fortunately however, the fact that (give or take random 
fluctuations) the curvatures are constant for systems in sta- 
tistical dynamical equilibrium (Eq. (1)) means that, for such 
systems, it is possible to get an estimate of these quantities by 
simply calculating them for different (pseudo) random realiza- 
tions of the same one particle phase space distribution. This 
in turn is sufficient to give us an idea of the values of the cur- 
vatures for fairly large spherical AT-body systems in statistical 
dynamical equilibrium. 

For systems consisting of up to A'^ = 1400, we integrated the 
full equations of motion and calculated the Ricci and scalar cur- 
vatures along the motion. This is done at intervals of a hundred 
starting for N = 100. The behaviour of the Ricci curvature 
time series during the evolution of such systems is described 
in El-Zant & Gurzadyan (1997). As in El-Zant & Gurzadyan 
(1997) we remove the contributions due to the closest encoun- 
ters by introducing a short range cutoff in the calculation of 
the potential energy and its derivatives. We choose a distance 
(including the softening) of 0.05 for this (i.e., 5% the virial 
radius). This is typically about a third of the minimum ra- 
dius of the neighbour spheres as calculated by the NBODY2 
algorithm. For systems consisting of up to A'' = 15000 parti- 
cles we calculate the Ricci and scalar curvatures for ten differ- 
ent realizations of the same equilibrium distribution and take 
the average. This is done starting at AT = 2000 at intervals 
of a thousand. To compare the results of the two approaches 
we take the average over the first ten outputs of the low N 
runs. This corresponds to one crossing time. Over this time 
interval the structure of the spheres are not expected to have 
changed much and, if indeed the motion is chaotic over small 
time-scales, these numbers also correspond to pseudorandom 
realisations of the same phase space distribution. The soft- 
ening of the potential (which also follows the Plummer law) 
was taken as 6 = 2/N which scales like the ratio of minimum 
to maximum impact parameters of standard relaxation theory 
(Binney & Tremaine 1987; Farouki & Salpeter 1994; Gierz & 
Heggie 1994). This is calculated as to exclude encounters (as- 
sumed independent) which lead to deflections beyond a certain 
maximal value. 

Except in the case of N = 100 when fluctuations are fairly 
large, causing it to be occasionally positive, it is found that 
the Ricci curvature (calculated while excluding the contribu- 
tion from neighbouring members as described above) is always 
negative. It is therefore easy to extract an average exponen- 
tial divergence time-scale through the method described in the 
previous section. The results are shown in Table 1. In the first 
column are the particle numbers. In the second column and 
fourth columns are the exponentiation time-scales (in cross- 
ing times) averaged over ten different values of the Ricci or 
Scalar curvatures as described above, while the third and fifth 
columns contain estimates of the RMS relative dispersion in 
the ten calculated exponentiation time-scales. This is obtained 
from 



where fe denotes the mean of the ten values. 

Two things are immediately clear from these results. The 
first is that the exponentiation time-scale is quite short — be- 



4 



iV 










100 


0.28 


0.78 


0.21 


0.06 


200 


0.34 


0.46 


0.23 


0.04 


300 


0.32 


0.13 


0.27 


0.01 


400 


0.28 


0.21 


0.29 


0.01 


500 


0.31 


0.17 


0.29 


0.01 


600 


0.28 


0.07 


0.29 


0.01 


700 


0.31 


0.12 


0.29 


0.01 


800 


0.31 


0.09 


0.30 


0.01 


900 


0.29 


0.07 


0.30 


0.01 


1000 


0.31 


0.09 


0.30 


0.01 


1100 


0.30 


0.07 


0.31 


0.01 


1200 


0.31 


0.05 


0.31 


0.01 


1300 


0.32 


0.05 


0.31 


0.01 


1400 


0.31 


0.05 


0.31 


0.01 


9nnn 


n 9Q 


n no 


U.Oi- 


n nn 

U.UU 


3000 


0.31 
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0.31 


0.00 


4000 


0.30 
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0.31 


0.00 


5000 


0.31 


0.00 


0.31 
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10000 


0.31 


0.00 


0.31 


0.00 
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0.31 


0.00 


0.31 


0.00 



Table 1. Variation with particle numbers of the exponentia- 
tion timc-scalcs Ter and Tes calculated from the values of the 
Ricci and scalar curvatures. In the third and and fifth columns 
the associated errors estimated with the aid of Eq. (8) are given 



Section 5 below). However, although its time averaged contri- 
bution is very small, it can have large positive values and, for 
small A'^, it occasionally causes the Ricci curvature to be posi- 
tive if no short range cutoff is introduced. When such a cutoff 
is introduced however, this term is always small compared to 
the third. It also becomes small compared to the third term 
when the particle numbers are large, thus supporting the pre- 
dictions of Gurzadyan & Savvidy (1984,1986) and Kandrup 
(1990a, 1990b) that as N increases the curvature is more likely 
to become negative. This can be seen from Table 2 where the 
values of the two terms are shown for single (pseudo) random 
realizations of Plummer spheres consisting of up to 45 000 un- 
softened particles (and no short range cutoff in the calculation 
of the potential and its derivatives) . 
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Table 2. Magnitudes of the second and third term in Eq. (1) 
and of the resulting Ricci curvature for large N . Values are ob- 
tained from a single pseudorandom realizations of the Plummer 
models 



ing a fraction of the crossing time. The second is that there 
appears to be no increase in the exponentiation time-scale as 
N increases. This is true of both the estimates using the Ricci 
and the scalar curvatures. Although the exponentiation time- 
scales predicted by the latter quantity slightly increase with 
when this number is limited to a few hundred particles, they 
quickly converge to the same value predicted by the Ricci cur- 
vature for larger A'^. 

Perplexing as these results may seem, their numerical ex- 
planation is actually relatively straightforward. First, since the 
Ricci curvature is mainly determined by the bulk properties 
of the system and not by fluctuations due to nearest neigh- 
bours, one expects its average over A'^ to remain constant with 
(large) increasing N, provided that these global properties (de- 
scribed by the one particle distribution function of the system) 
remain unchanged. This means that the exponentiation time- 
scale, which depends on that average (cf. Eq. (6)), remains 
constant for large enough A'^. The order of magnitude of the 
exponentiation time-scale can be understood as follows. The 
last term on the right hand side of Eq. (1) is very small when 
the softening is the situation we are interested in here. The 
first term is also small if the forces on the particles are not 
aligned with their velocities (a very improbable situation in 
an equilibrium configuration). We are therefore left with the 
second and third terms in the expression for ru. For systems 
near virial equilibrium, and in the absence of ordered motion, 
the second term is dominated by contributions from close en- 
counters. It is highly fluctuating and averages to zero when 
the softening is small (this statement however does not ap- 
pear to hold when ordered motion is present as we shall see in 



The above explains why the exponentiation time-scales de- 
scribed by the Ricci and Scalar curvatures are similar since it 
is the third term that appears in the formula of the scalar cur- 
vature (e.g., Gurzadyan & Savvidy 1986). A simple argument, 
due to Kandrup (1989), shows why these time-scales should 
not vary much with particle number. Adopted for use in con- 
junction with the Ricci curvature, it goes as follows. For large 
AT the Ricci curvature dominated by the last term becomes 

where F denotes the average total RMS force per unit mass 
acting on a test particle. Using Eq. (6), this implies an expo- 
nentiation time-scale of Te ~ u/F, where v is the RMS speed. 
This time-scale is of course of the order of a dynamical time. 

4. Behaviour of the curvatures as the softening param- 
eter is increased 

In the case of A'^ = 1000, we have calculated the Ricci and 
scalar curvatures for a range of softening radii starting from 
10^''' units (i.e, one thousandth of the virial radius) to 4 x lO^'^ 
units. In Fig 1 we plot the Ricci curvature as a function of the 
crossing time for some representative values of the softening 
parameter. These plots show that as this parameter increases 
so does the Ricci curvature, eventually it becomes positive just 
as in the case of the flat systems (El-Zant 1996a). 

In Table 3 are shown the exponentiation time-scales cal- 
culated from the values of the Ricci and scalar curvatures av- 
eraged over the first crossing time for the various simulations 
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Table 3. Variation of the exponentiation time-scales with the 
softening parameter tp. Other symbols are as in Table 1. The 

Ricci and scalar curvatures arc positive for the last two values. 
For this reason no exponentiation time-scale is defined 



that were run. Since the exponentiation time-scale has an in- 
verse square root dependence on the Ricci curvature, it is relar 
tively insensitive to this quantity unless it is close to zero. The 
range of softening parameters in which the transition occurs 
from negative values of — corresponding to an exponentia- 
tion time-scale of one dynamical time — and positive ones is 
rather small. Therefore, although in principle arbitrarily large 
exponentiation time-scales can be obtained by fine tuning the 
softening parameter, in practice a softened system appears to 
have an exponentiation time-scale of the order of a dynamical 
time or none at all. 

Prom Table 3 one can also see that, for large softening, the 
exponentiation time-scales associated with the Ricci curvature 
are rather large compared to those derived from the average 
value of the scalar curvature. This is because the second term 
in Eq. 1, which is not present in the formula for the scalar cur- 
vature, is always positive when significant softening is present. 
In turn this is because this term has contributions which con- 
tain the second derivatives in the form 

(10) 

In a spherically symmetric system with positive density de- 
creasing towards the outside, this expression is positive. It can 
be fairly large when the softening is large. However it is still 
highly fluctuating — hence the large errors in the exponentiar 
tion time-scales. This is of course in addition to the last term 
in (1) which is also present in the formula for the scalar curva- 
ture and is always positive for systems with positive density. 

A completely analogous situation transpires when one in- 
creases N with fixed softening and no short range cutoff. Inde- 



pendent of the softening length or its precise functional form, 
the Ricci curvature is positive and its average over N remains 
constant as one varies the particle number. This is easy to un- 
derstand since particles which will contribute significantly to 
the local density at a given point are those within a distance 
of the order of magnitude of the softening length h. For high N 
the number of particles in such spheres increases as Nb^ while 
the contribution to the density from each of the particles goes 
down as 1/6^. 

At first sight, the above results may be taken to mean that 
the negative Ricci curvature is somehow caused by the singu- 
larity in the Newtonian potential or due to the contribution of 
nearby neighbours — and again the predicted instability may 
be a mathematical artifact. A closer look however reveals that 
this is not so. Near the transition from negative to positive cur- 
vature, the assumptions justifying the validity of the negativ- 
ity of the Ricci curvature as an indicator of average instability 
break down (since in this case many of the two dimensional 
curvatures will already be positive and the variation in their 
absolute values may be very large) . The positivity of the Ricci 
curvature in this case only means that the relative motion of 
nearby particles is regular since the (large) force between them 
is approximately constant with distance. These contributions 
dominate the value of the Ricci curvature. A similar situation 
for example occurs if one encloses the system in an "elastic 
sphere" where particles near the boundary are subjected to a 
harmonic potential. If the spring constant is large, the curva- 
ture is dominated by the resulting positive contributions from 
these particles (this experiment was actually conducted by the 
author). Obviously the nature of the gravitational dynamics in 
this case is not radically modified. 

The transition to positive curvature therefore should not 
be viewed as indicating a sharp switch from a chaotic to an 
integrable system but from one where the majority of trajec- 
tories were highly unstable to one where their instability is 
somewhat less pronounced (how pronounced can only be de- 
termined by calculating the two dimensional curvatures and 
integrating the full set of linearised equations). Only in the 
true continuum limit (infinite N and fixed softening) is full 
integrability recovered (it is interesting to note that the equiv- 
alence of the continuum limit to the large- A'^ limit can only be 
proven for twice diffcrcntiablc potentials which are bounded ev- 
erywhere: Braun & Hepp 1977; Spohn 1980). Nevertheless, the 
fact that the curvature is affected by softening does suggest 
that the origin in the instability of the trajectories of gravi- 
tational systems is related to their discrete nature — this is 
in line with the fact that smoothing out the force field elimi- 
nates this instability. The eff'ect of softening on the exponenti- 
ation time-scale was studied directly by Suto (1991). He found 
that the exponentiation time-scale was related to the softening 
by Te ~ 206^^^rc (where h is in units of average interparticle 
distance). Thus, while the softening docs increase the relax- 
ation time-scale, it does it in a rather moderate manner. On 
the other hand, only including the contribution due to nearby 
neighbours in the calculation of the curvature yields an expo- 
nentiation time-scale that increases as iV^/^Tc (Gurzadyan & 
Savvidy 1986). We therefore conclude that while the trajec- 
tory instability in gravitational systems appears to be related 
to their discrete nature, what gives rise to this property should 
be the contribution to discreteness noise from the whole sys- 
tem and not just from neighbouring particles. We discuss this 
further in the concluding section. 
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5. The effect of rotation 

In El-Zant (1996a) it was found that self gravitating sheets 
starting from states of sohd rotation initially had a larger av- 
erage Ricci curvature than ones starting from random initial 
conditions in velocity space. This property is important since 
the fact that many flattened disk gala^xies consist predomi- 
nantly of stars moving on nearly circular trajectories near the 
disk plane suggests that for such systems phase space diffusion 
time-scales must be relatively long. In this section we study in 
more detail the change in the Ricci curvature as the energy in 
rotational motion is increased. 

Constructing stable equilibrium spherical self gravitating 
systems with a wide range of ratios of rotational to random 
motion and having the same density distribution is not a triv- 
ial task (e.g., Palmer 1994), we therefore stick to the simple 
situation of static averages. The systems chosen here are ho- 
mogeneous and — in the absence of rotation — have isotropic 
velocity distribution which does not vary with radius. Solid 
body rotation is then added to the random motion before all 
velocities are rcscaled so as to have a total energy of -0.25 in ac- 
cordance with the system of units discussed above (Section 3). 
As we have done before, we calculate the Ricci curvature for ten 
different pseudorandom realizations of the same distribution, 
calculate the average, and estimate the RMS error (since the 
scalar curvature docs not explicitly depend on the velocities it 
remains unchanged when the velocity distribution is changed). 



Trot/T 




er 


0.00 


0.30 


0.07 


0.02 


0.43 


0.08 


0.15 


0.46 


0.08 


0.33 


0.51 


0.08 


0.55 


0.61 


0.08 


0.70 


0.71 


0.07 


0.81 


0.84 


0.07 


0.90 


1.02 


0.06 


0.95 


1.18 


0.05 


0.97 


1.26 


0.04 


0.98 


1.32 


0.04 


0.99 


1.39 


0.03 


0.995 


1.42 


0.03 



Table 4. Variation of the exponentiation time rer, averaged 
over ten (pseudo)random realizations of N — 1000, as the frac- 
tion of kinetic energy of rotational motion Trot/T is increased. 
er is an error estimate obtained by using Eq. (8) 

Table 4 shows the variation of the corresponding exponen- 
tiation time-scales as the energy in rotational motion is in- 
creased. As is clear from these results, the exponentiation time- 
scales are significantly increased when the rotational motion is 
increased. Looking at the different terms on the right hand 
side of Equation (1), it is easy to see that only the first two 
are directly dependent on the velocities and may therefore be 
affected by the reordering of random motion into rotational 
motion. Of these two terms the first — as discussed above — 
is small for most equilibrium velocity distributions. It is even 
smaller for rotating systems since the sum of the scalar product 



of the particles' velocities and the forces acting on them which 
this term represents is near zero. The second term consists of 
two parts. One involves the quantities in (4) and represents the 
second derivatives of the potential with respect to coordinates 
of the same particle multiplied by the velocities of that parti- 
cle. This term is also relatively small since again the gradient 
of the force in a corresponding smoothed out system is nor- 
mal to the velocities (albeit not as small as the first term since 
the components of the derivatives of the force have larger fluc- 
tuations than the components of the force) . The other part of 
the second term (involving the quantities in (3)) consists of the 
derivatives of the force at one particle's position with respect to 
another's projected on the velocities of the two particles — in 
other words it measures correlations in velocities and positions 
between the trajectories of particles in the systems. This term 
is small and fluctuating when the velocity field is random but 
is much larger and has positive sign when the kinetic energy is 
in the form of ordered rotational motion. 

The fact that the exponentiation tinic-scalo is significantly 
larger when ordered motion is present confirms that, except 
when effects arising from the non-compactness of the phase 
space are important (e.g., escape of particles), higher thermo- 
dynamic entropies appear to be related to higher values for the 
dynamical entropy. Thus confirming that A'^-body systems will, 
in general, evolve towards higher dynamical entropy states as 
was found to be the case in El-Zant & Gurzadyan (1997). 

If directly interpreted to mean mixing on an exponential 
rate as is expected in the standard case of gravitational sys- 
tems with negative two dimensional curvatures, the derived 
time-scales predict that a rotating system will evolve on a time- 
scale of about 45 crossing times (El-Zant 1996a). This is still 
compatible with the age of average disk galaxies at about 10 
kpc, and it is possible that for realistic density and velocity 
distributions the predicted evolution time-scale may be still 
larger. In addition, enough two dimensional curvatures may 
be positive so as to restrict the motion. Therefore, as in the 
case of softened systems, the variation of the Ricci curvature 
exponentiation times with rotation should only be interpreted 
as representing a trend. However, as we shall see in the next 
section, it appears that the relation between the exponentia- 
tion time-scales and the macroscopic evolution of gravitational 
systems may not be so direct. 

6. Conclusions and possible interpretation 

In this paper we continued our investigation of the behaviour of 

the Ricci and scalar curvatures of the configuration manifolds 
of A''-body gravitational systems. These relate the geometry of 
the phase space to the stability properties of trajectories on 
it. It was found that, for spherical systems with isotropic ve- 
locities, the inferred exponentiation time-scale is rather short 
(less than a crossing time) and did not depend on the parti- 
cle number when the softening length decreased as 1/N and 
a short range cutoff in the potential was introduced. The ex- 
ponentiation time-scale however was found to be affected by 
the presence of ordered rotational motion or when the softening 
was increased (while keeping the particle number fixed). In the 
first case it was found to increase significantly while in the sec- 
ond case it could even become undefined because the curvature 
became positive. A similar process takes place if the softening 
radius is fixed and the number of particles is increased (if no 
short range cutoff is introduced). 
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The exponentiation time-scale being so short and not vary- 
ing with particle number means that it is difficult to uncover 
its significance. The main problem is that these properties ap- 
parently contradict the intuitive idea that particles in large 
systems should move essentially unperturbed in the mean field 
potential. In spherical potentials this happens to mean that 
they all lie on regular trajectories. It then should follow that 
the divergence between nearby trajectories is, on average, lin- 
ear and not exponential. 

One may like to relate the exponential instability to the 
process of achieving dynamical equilibrium, the time-scale of 
which does not vary with particle number, (e.g., Kandrup 1989; 
El-Zant 1996c). While this could be the case, one expects 
that even completely smooth spherical systems achieve such 
an equilibrium. The exponential instability on the other hand 
appears to be inherently related to the discreteness of iV-body 
gravitational systems. At the same time this does not imply 
that it is mainly a result of close encounters as was explained 
near the end of Section 4. In the light of that discussion, and 
looking at the relative importance of the terms of the formula 
for the Ricci curvature (see Section 3), one comes to the con- 
clusion that, for spherical systems with isotropic velocity dis- 
tribution and in virial equilibrium, the negative curvature is 
related to the fact that A'^-body systems have a largo mean 
field force (because the interaction is long range) and at the 
same time have locally peaked density distribution (because of 
their being composed of particles). Thus we may expect that 
the cause of the exponential instability is the discreteness ef- 
fects due to the long range full A-body interaction. 

To see how long range gravitational interaction can trig- 
ger chaotic behaviour, we follow Chirikov (1979) and divide 
the Hamiltonian of the system under consideration into an 
unperturbed part Ho and an non-integrable time dependent 
perturbation V. In terms of action angle variables this reads 

H{J,@,t) = Ho{J) + ViJ,&,t). (11) 

For our purposes Hq will be the smooth spherically symmetric 
potential and V would the perturbation arising from discrete- 
ness effects: V = Vnbody — Vb. If we assume that V is periodic 
in time with phase r = Q,t + To, there will be a whole set of 
resonance conditions given by 

miUJi{3i)+nn = Q, (12) 

where uji {i = 1,3) are the natural frequencies of oscillation 
of the action variables of the motion in the background po- 
tential Vi). Around the neighbourhood of these resonances a 
stochastic layer on which chaotic motion can take place forms 
(e.g., Lichtenberg & Lieberman 1983 (LL))- If ^ >> u)i how- 
ever, the resonance conditions are satisfied only for very large 
m/n and no lower order resonances occur. In this case the ef- 
fect of these resonances is limited because the stochastic layer 
around them is small (LL; Moiss 1987). The vast majority of 
non-resonant trajectories will therefore remain stable. In the 
six dimensional phase space this will mean the perturbation 
causing trajectories to move between KAM tori. In that case 
the original idea of Chandrasekhar (1942) is justified: he con- 
sidered discreteness effects duo to short lived encounters with 
nearest neighbours where resonance effects are negligible and 
where there is a clear separation of time-scales justifying the 
assumption of independence. However it is now thought that 



weak distant encounters dominate two body interactions in N- 
body systems and that distant encounters are more important 
for larger systems (BT). Thus, although the strength of per- 
turbations due to discreteness decreases as ~ the den- 
sity of resonances (per action per particle) increases as ~ AT 
and is increasingly dominated by more effective terms. This 
might explain the observed persistence of chaotic behaviour 
for large N. For, according to LL, this issue reduces to "the 
question of whether the density of important resonances, as 
projected at a single action, increases faster than the width of 
the resonances decreases. If this happens then we would expect 
resonance overlap and strongly chaotic motion to occur for N 
degrees of freedom as A ^ oo" . 

It is important to note here that although the main cause 
of the chaotic behaviour may be distant interactions between 
particles in a A?^-body system, this does not mean that the 
exponential instability will not be affected by short range en- 
counters. For, even if one considers these as additional high 
frequency noise added to the system, according to Pfenniger 
(1986) such perturbations completely change a chaotic sys- 
tem's trajectory on relatively short time-scales. One therefore 
expects high froquroncy discreteness noise to be an additional 
source of instability which will affect the exponentiation time- 
scale. This would explain why, even if the instability is mainly 
caused by distant encounters, the exponentiation time-scale in- 
creases when the force softening is increased this suppressing 
the high frequency noise. 

It is clear that as N increases, quantities such as the en- 
ergy and angular momenta of individual particles will be better 
conserved. This however does not imply that the decorrelation 
time-scale defined by the exponential divergence should be- 
come smaller. For example a perturbed pendulum can display 
highly chaotic oscillations which decorrelate very fast while 
changes in the amplitudes of these oscillations are much slower. 
In this type of situation one can average over the fast phase 
variable which may (because of the short decorrelation time) 
be considered as random. The evolution of the actions can then 
be regarded as a diffusion process. In fact this approximation is 
only valid when the dynamics is strongly chaotic so that chaos 
occurs for the vast majority of initial conditions (Chirikov 1979; 
LL; Shlesinger et al. 1993). This is of course what appears to 
be the case for gravitational A-body systems. 

In the case of gravitational systems, it is actually not that 
surprising that the exponentiation time-scales are very differ- 
ent from the diffusion time-scales of the action variables. For 
in the standard case of a system with negative two dimensional 
curvatures, the exponential instability will guaranty that the 
system visit all regions of the available phase space (Anosov 
1967). In the case of a A-body system this will be the whole 
subspace defined by the conservation of total energy and mo- 
mentum. Since in a Hamiltonian system the phase density is 
conserved along the motion, this will imply that the phase 
space (A^ particle) distribution function will become constant 
when averaged over progressively smaller volumes in that space 
— that is the system becomes more and more likely to be found 
in any of its microscopic states. For an open gravitational sys- 
tem this cleaxly cannot be the case since, instead of evolving to- 
wards a definite thermodynamics equilibrium, this type of sys- 
tem contirmally evolves towards more and more concentrated 
configurations when it divides into a contracting core section 
surrounded by an expanding halo. In this type of evolution, the 
Poincare recurrence theorem is not valid and the system need 
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not, even in principle, return to less inhomogeneous states. In- 
stead it continually moves into new areas of the phase space 
characterised by larger entropy. In this case, if the evolution 
time was directly related to the exponentiation time-scale, then 
there would be no chance for the distribution function to be 
constant on any region of the phase space. Therefore, no type of 
equilibrium would be possible and gravitational systems would 
disintegrate in a few dynamical times! What prevents this from 
happening of course is that some states are long lived because 
they are stable dynamical equilibria. 

Indeed, it was found that, for closed systems which did 
have a definite thermodynamic equilibrium state, this state 
was reached on a time-scale comparable to the exponentiation 
time if no dynamical equilibrium existed between the initial 
state and the final equilibrium (this happened even when the 
virial ratio remained nearly constant during the evolution thus 
ruling out violent relaxation). On the other hand, if there ex- 
isted intermediate dynamical equilibrium states between the 
final equilibrium and the initial configuration, the relaxation 
to the final equilibrium state took much longer (El-Zant 1997) . 
It may be therefore that the exponential instability acts as to 
smooth out the distribution function on the subspaces where 
these equilibria are defined while at the same time driving the 
evolution towards higher entropy states. This means that al- 
though the instability is locally normal to the phase space mo- 
tion, the global stretch and fold of that space is such that the 
decorrelation is faster along the motion — justifying the diffu- 
sion description for the evolution of the action variables. For 
larger N, the dimension of the phase space grows and there are 
more states (corresponding to a given dynamical equilibrium) 
to be covered while the exponentiation time-scale is constant. 
Thus, one expects the diffusion rate of the action variables to 
be larger as N increases. 

More precisely, it can be shown using symbolic dynamics 
(Alexseev & Yakobson 1981; McCauley 1993; Chirikov 1994) 
that for hyperbolic systems (a class which includes systems 
with negative two dimensional curvatures), with well defined 
final states, the time-scale for a system to achieve such a state 
with coarse grained volume resolution is 

Tr - — (13) 

where h ~ (\/—2ru) is the Kolmogorov entropy in the CEise 
when Tu is dominated by the largest two dimensional curva- 
tures. (The above relation is exact and follows from equation 
(14.3) of Alexseev & Yakobson when h is constant along the 

Ai'-body trajectory, otherwise it is valid in an averaged asymp- 
totic sense). Fixing the linear resolution d and using Eq 6 one 
finds 

Tr = 2V3Hl/d)W^Te ~ ■ (14) 

Therefore, in systems for which definite statistical equilibria 
exist, the time-scale to achieve such equilibria docs indeed in 
general increase with A'^ and appears to scale as the inverse 
of the strength of the perturbation due to discreteness. This 
is the case even though the exponetiation (or Liapunov) time- 
scale as derived from the value of the Ricci curvature (or as 
directly inferred from the divergence of A''-body trajectories) 
does not vary much with N. This is essentially because, unlike 
the case of low dimensional systems, the Kolmogorov entropy 
is dominated by the maximal exponent and volume resolution 



becomes very different from linear resolution. Thus, for large- A^ 
systems, the Liapunov time-scale does not necessarilly coincide 
with the mixing time-scale. A more detailed discussion of the 
mathematical origin of this important effect, which in retro- 
spect has been a major cause of the confusion regarding the 
meaning of the Liapunov time-scale for A^-body systems, will 
be given elsewhere. 

Because of the exponential instability, the (6-d) phase space 
trajectories of particles quickly decorrelate so that the proba- 
bility of finding a particle in a given state is independent of the 
state of the other particles. This being so even if the motion is 
affected by the discrete nature of the system. For times longer 
than the decorrelation time Te the A^ particle distribution func- 
tion can then be given in terms of the individual one particle 
distribution function by 

= /(xi, xi, i)/(x2, X2, i).../(x]v, xjv, t). (15) 

This is of course the requirement that a system be completely 
described by the coUisionless Boltzmann equation (CBE). 
Moreover, since steady state solutions of the CBE will only 
depend on the action variables (Jean's theorem), these solu- 
tions will be equivalent to those produced by the mean field 
dynamics — for times long compared to the exponentiation 
time-scale but short compared to the diffusion rate of the ac- 
tion variables. In this context the continuum (coUisionless) ap- 
proximation may then be valid in an averaged sense: the exact 
trajectories in smoothed background potentials would not be 
valid but the time averaged orbits would be correct as long as 
the time-scales considered are small compared to the diflfusion 
times of the action variables. 

Although for many purposes the situation described above 
is similar to that of standard galaxy dynamics, it diff'ers in 
one important respect: the underlying motion is intrinsically 
chaotic and cannot be expressed as a linear superposition of reg- 
ular motion and independent binary encounters. For this reason 
the diffusion time of the action variables need not be accurately 
represented by the two body relaxation time. Moreover, it is 
now well known that chaotic trajectories can repond to exter- 
nal perturbations in a manner that is different from that of 
regular trajectories (Pfenniger 1986; Kandrup 1994; Merritt &: 
Valluri 1996; El-Zant 1996c). Therefore, the response of the 
trajectories to additional perturbations (e.g., high frequency 
discreteness noise or global asymmetries in the background po- 
tential) is not necessarilly identical to that of the trajectories 
in the smoothed out potential — which for spherical systems 
happen to be all regular. 

To sum up. In this section it has been suggested that reso- 
nances between the orbital frequencies of particles and forcing 
caused by the dicreet nature of the global potential give rise 
to chaotic trajectories which decorrelate over a time-scale of 
the order of a dynamical time. Nevertheless some quantities 
may decorrelate on the much larger time-scale. In the case of 
systems of having a well defined final equilibrium state the 
maximum such time-scale is that needed to reach a statistical 
equilibrium when such a state exists which scales as Af^^^. 

The above ideas can be tested in more than one way. As 
mentioned in the introduction, it is well documented that AT- 
body trajectories decorrelate over a dynamical time or less. It 
may be useful however to also examine the stability of individ- 
ual particle trajectories in A'-body simulations — preferably 
by methods that do not require linearisation of the equations 
of motion (e.g., Laskar's frequency analysis). To find out the 



9 



dominant range of encounters causing the exponential diver- 
gence in A'^-body systems, it may be possible to calculate the 
frequency spectrum of perturbations a A^-body particle is sub- 
jected to and examine the effect of its different regions on in- 
dividual trajectories in smoothed out potentials. Alternatively, 
one may want to bin the contributions to the discreteness noise 
F — Fsmooth — Fnbody actiug along a particle's trajectory into 
impact parameter ranges and again examine their effect on the 
stability of test particles' trajectories in the smoothed out den- 
sity distribution. The time-scales over which the physical char- 
acteristics of a given system changes can be studied directly by 
examining the relaxation of the particles' integrals of motion 
in the smoothed out potential. This can be done either by com- 
puting the diffusion coefficients of these variables for particles 
in A'^-body simulations or by studying the macroscopic relax- 
ation of numerically simulated systems. The latter approach is 
most effective for closed systems where definite thermal equi- 
libria are well defined (El-Zant 1996b). One would then look 
at the time-scale for attaining isothermal equilibrium and the 
time-scale of relaxation of initial anisotropic velocity distribu- 
tions (El-Zant 1997) and how these time-scales vary with A'^ 
(El-Zant & Goodwin 1997). In fact, the existence of a well de- 
fined final state that ceases to evolve also means that all the 
aforementioned tests are easier to conduct for closed systems. 
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Figure caption 



Fig. 1. Evolution of the Ricci curvature for N = 1000 systems with 
diflferent values for the softening parameter b which takes the values 
of (from top to bottom) 4.0 X 10-^, 1.6 X 10-^, 2.8 x 10"^ and 
4 X 10-2 
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